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boundary  or  transition  region  between  the  high  and  low  resolution  solutions. 
For  any  problem  of  practical  interest  there  is  no  "correct"  wav  of  handling 
the  boundary  problems  and  no  way  of  avoiding  some  boundary  generated  noise. 
, ^t'onswtwe-fttiy  u pragmatic  approach  is  proposed  here  which  overspecifies  the 
artificial  lateral  boundaries  of  a fine-mesh,  limited-area,  nested  region  and 
depends  upon  a high-order  phase  and  amplitude  restoring  interpolation  along 
with  a high-order,  but  weak,  smoothing  operator  to  control  noise  and  ensure 
verisimilitude.  ) 

The  basic  problem  as  well  as  the  method  proposed  is  outlined  in  some 
detail  in  the  Introduction.  High  order  amplitude  restoring  and  phase  and 
amplitude  restoring  interpolation  operators  are  defined  and  constructed  in 
Section  2,  In  Section  3,  both  static  and  dynamic  comparisons  are  made  be- 
tween the  two  types  of  interpolation  operators.  In  the  static  tests,  the  phase 
and  amplitude  restoring  interpolation  is  shown  to  be  superior  in  minimizing 
root  mean  square  differences  between  interpolated  values  and  analytic  values 
of  simple  Fourier  functions. 

The  dynamic  tests  are  carried  out  in  the  framework  of  a multilevel 
| primitive  equation  model  and  show  unequivocally  the  importance  of  correct 
phase  information  for  the  lateral  boundaries  of  the  fine-mesh  region.  The 
results  obtained  indicate  that  the  simple  procedure  proposed  here  should  yield 
satisfactory  results  in  operational  fine -mesh,  limited-area  models. 
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The  Treatment  of  Lateral  Boundary  Conditions  in 
Limited -Area  Models:  A Pragmatic  Approach 

I.  INTRODUCTION 

Routine  prediction  of  atmospheric  behavior,  particularly  the  prediction  of  the 
gross  features  of  the  large-scale  circulation,  is  commonly  based  upon  the  numeri- 
cal solution  of  some  mathematical  model  which  is  designed  to  incorporate  the  phy- 
sics for  those  atmospheric  features  of  immediate  interest.  The  typical  mathemati- 
cal model  consists  of  a set  of  nonlinear  partial  differential  equations  expressing  the 
rate  of  change  of  momentum  and  heat  (the  prognostic  equations)  and  an  appropriate 
set  of  diagnostic  equations  expressing  relationships  among  the  prognostic  variables 
and  other  variables  which  are  required  in  the  numerical  solution  of  the  prognostic 
equations. 

The  system  is  solved  numerically  as  an  initial  value  problem  over  a network  of 
grid  points  in  space.  The  distance  between  neighboring  points  in  this  network  con- 
trols the  detail  that  can  be  resolved  in  the  predicted  fields.  However,  though 
arbitrary  to  some  extent,  this  distance  is  more  or  less  determined  by  outside  factors 
such  as  the  scale  size  of  the  phenomena  of  interest,  but  especially  by  the  available 
computer  capacity.  The  required  computer  capacity  is  particularly  sensitive  to  the 
grid  distance  since  a halving  of  the  mesh  size  in  a three-space  dimensional  model 
increases  the  computer  time  by  more  than  a factor  of  10,  assuming,  optimistically, 
that  the  computer  has  sufficient  unused  internal  (high  speed)  storage  capacity  to 
avoid  a concomitant  increase  in  the  book-keeping  operations. 

(Received  for  publication  18  April  1977) 
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It  is  often  desirable  to  increase  the  forecast  resolution  over  some  limited  area 
of  the  globe  where  the  observational  density  or  the  nature  of  the  problem  warrants 
increased  resolution,  but  where  economic  factors  make  it  impractical  to  increase 
the  resolution  uniformly  over  the  whole  globe.  in  tact,  even  if  the  machine  capa- 
city were  available,  it  might  be  desirable  in  the  interest  of  speed  and  economy  to 
use  a coarse-mesh  resolution  in  regions  where  the  fields  are  smooth  and  changing 
slowly  and  to  use  a fine-mesh  resolution  only  in  regions  where  the  fields  are 
changing  rapidly  in  time  and  space.  In  such  cases  the  problem  no  longer  assumes 
the  character  of  a pure  initial  value  problem,  but,  depen  ling  upon  the  numerical 
proce  lure,  assumes  more  or  less  the  character  of  an  initial  boundary -value  prob- 
lem. There  are  a variety  of  ways  in  which  such  a problem  can  he  handled. 

(It  Perhaps  the  simplest  is  to  choose  the  'ine-ntesh  domain  large  enough  so  that 
during  the  course  of  the  forecast  the  lack  of  knowledge  of  the  boundary  values  can- 
not propagate  into  an  interior  region  of  interest  and  affect  the  value  of  the  forecast 
in  this  interior  region.  This  approach  might  be  necessary  il  there  were  no  informa- 
tion outside  the  boundary  of  the  fine-mesh  domain.  Boundary  errors  are  simply 
accepted  and  are  kept  under  control  by  strong  smoothing  near  the  boundaries. 
Inasmuch  as  this  approach  requires  a large  fine-mesh  region  it  is  computationally 
costly  and  in  spite  of  its  simplicity  it  is  not  generally  used  if  there  are  reasonable 
alternatives. 

(2)  In  *he  variable  mesh  approach  a single  solution  is  obtained  for  the  'global' 
domain,  but  the  mesh  density  is  allowed  to  vary  so  that  regions  of  interest  that 
warrant  increased  resolution  are  covered  bv  a fine-scale  mesh,  whereas  the  re- 
mainder of  the  domain  is  covered  by  a eoarse-seale  mesh.  The  principal  advantage 
of  this  approach  is  that  there  are  no  artificial  lateral  boundaries  around  the  fine- 
mesh  domain  for  which  separate  boundary  conditions  have  to  be  specified.  Further- 
more, this  approach  permits  a two-way  interaction  between  the  flow  fields  defined 
by  the  fine-and  coarse-mesh  regions. 

2 

This  approach  was  used  by  Birchfield”  in  a moving  coordinate  system  for  a 
hurricane  problem.  Moretti^  suggests  a stretching  transformation  of  the  coor- 
dinates which  lias  the  virtue  of  permitting  smooth  transitions  between  the  regions  of 

different  grid  lensitv.  The  variable  mesh  approach  has  been  applied  by  Harrison 

4 .5.6 

and  Klsberrv  to  a two-dimensional  multilevel  model  and  by  Koss,  ‘ Ookochi.  and 

Phillips  and  Shukla1  to  the  shallow -water  equations.  The  variable  mesh  method  was 

adopted  by  Mathur^  in  a four-level  primitive  equation  hurricane  model,  while  Price 

and  Mac Pherson  ' discussed  the  use  of  cubic  polynomial  splines  to  fit  the  spatial 

variations  of  the  dependent  variables. 

(Because  of  the  large  number  of  references  cited  above,  they  will  not  he  listed  here. 
See  Reference  Page  No.  67,  for  References  1 through  hi. 
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In  spite  of  the  advantages  of  this  approach,  as  indicated  above,  the  results 

have  been  generally  disappointing  except  in  eases  where  the  large-scale  domain 

fields  tend  to  be  smooth  "ith  weak  gradients  (such  as  in  the  hurricane  problem),  or 

where  smoothness  is  ensured  by  the  numerical  technique  (such  as  with  the  use  of 
n 

bicubic  splines.  The  inherent  difficulty  with  the  variable  mesh  approach  arises 
from  the  fact  that  the  phase  speed  of  the  waves  is  a function  of  the  truncation  error, 
which  in  turn  lepends  upon  the  grid  spacing  as  well  as  upon  the  dominant  scale 
sizes  of  the  field  variables.  Waves  in  the  eoarse-mesh  domain  will  change  phase 
speed  when  passing  into  the  fine-mesh  region  and  will  interact,  producing  inter- 
ference with  the  part  of  the  wave  that  has  remained  in  the  eoarse-mesh  domain. 

As  Browning  et  al'*'  point  out,  this  interference  phenomenon  is  present  regardless 
of  the  difference  method  used. 

It  seems  intuitively  apparent  that  abrupt  changes  in  grid  density  would  produce 
more  serious  interactions  than  a gradual,  nhased  transition.  Nevertheless,  in  view 
of  the  increased  complexity  in  programming  necessitated  by  this  approach  and  in 

view  of  the  dubious  advantages  (inasmuch  as  changes  in  grid  density  induce  effects 

11 

similar  to  those  produced  bv  artificial  lateral  boundaries,  it  appears  worthwhile 
to  investigate  a third  general  method. 

(3)  The  third  approach  entails  separate  eoarse-mesh  and  fine-mesh  calculations, 
with  the  fine-mesh  region  nested  within  the  eoarse-mesh  domain.  The  eoarse- 
mesh  solution  is  obtained  for  a global’  domain  without  artificial  lateral  boundaries 
and  serves  to  supply  the  necessary  boundary  information  for  a solution  with  a fine- 
mesh  grid  over  a limited  area  of  the  globe,  The  principal  difficulty  with  this  ap- 
proach concerns  the  choice  of  proper  boundary  conditions  for  the  particular  svste  n 
of  equations.  This  problem,  while  serious,  appears  amenable  to  treatment  ami  the 
nested-domain  approach  seems  fundamentally  simpler  for  routine  application  than 
the  variable-grid  approach. 

The  choice  of  proper  boundary  conditions  for  the  nested  region  has  been  exten- 

12-14 

sively  examined  for  simple,  linearized  hyperbolic  systems.  The  boundary  a 

conditions  should  provide  a well -posed  problem.  That  is,  the  solution  should  de- 
pend continuously  upon  the  boundary  data  so  that  a small  change  in  a boundary  value 

10.  Browning,  Cl.  , Kreiss,  11. -O. , and  Oliger,  J.  (1073)  Mesh  refinement.  Math. 

Comp.  , 27:20-39. 

W«A 

11.  Kreiss,  H. -O.  , and  Oliger,  ,1.  (11173)  Methods  for  the  Approximate  Solution  of 

Time  Dependent  Problems,  CJARP  Pub.  Scr.  NIo.  TTTj  103  pp. 

12.  Kreiss,  H.  -O.  (1070)  Initial  boundary  value  problems  for  hyperbolic  systems. 

Comm,  Pure  and  Appl.  Math.  , 23:277  -208. 

13.  Gustafsson,  B.  , Kreiss,  H. -O.  , and  Sundstrom,  A.  ( 1972)  Stability  theory  of 

difference  approximations  for  mixed  initial  boundary  value  problems,  11. 

Math.  Comp.,  26:640 -686. 

1—  vW 

14.  Klvius,  T.  , and  Sundstrom,  A.  ( 10731  Computationally  efficient  schemes  and 

boundary  conditions  for  a fine-mesh  barotropic  model  based  on  the  shallow- 


should  have  only  an  appropriately  small  effect  on  the  solution.  However,  the  non- 
viscous  baroelinic  primitive  equations  appef  to  be  ill-posed  for  any  specification 

15 

of  the  lateral  boundary  conditions.  The  difficulty,  which  is  due  at  least  in  part 

to  the  hydrostatic  approximation,  6 can  be  alleviated,  if  not  overcome,  by  the  ex- 

11  I r> 

pedient  of  an  artificial  viscosity.  ' Additional  boundary  difficulties  are  intro- 
duced in  solving  the  finite-difference  equations,  especially,  when  for  purposes  of 
reducing  the  truncation  error,  the  finite -difference  equations  are  of  higher  order 
than  the  original  differential  equations,  thus  requiring  computational  boundary  con- 
ditions1 and  consequently  giving  rise  to  computational  as  well  as  physical  solu- 

, . It! 
tions. 

Although  unable  to  develop  a general  theory  on  the  existence  and  uniqueness  of 

15  14 

solutions  to  nonlinear  hyperbolic  systems,  Sundstrom  and  Flvius  and  Sundstrom 
discuss  sufficient  conditions  for  the  well -posedness  of  the  pure  initial  value  prob- 
lem for  the  shallow -water  equations.  For  certain  a -posteriori  boundary  conditions, 
the  well  -posedness  also  applies  to  the  limited  area  problem  (that  is,  to  the  initial- 
boundary value  problem).  Although  it  might  appear  reasonable  to  expect,  it  is  not 
possible  to  state  with  confidence  that  for  a given  initial  state  the  solution  wi'c  re- 
main continuous  and  differentiable. 

In  the  case  of  the  finite-difference  equation  for  the  limited-area  problem, 

there  remains  even  greater  uncertainty.  The  stability  of  finite-difference  solutions 

jq  20 

has  been  shown  for  simplified  linear  systems  bv  the  use  of  the  energy  method.  ’ “ 

14 

For  the  general  nonlinear  limited-area  problem  Elvius  and  Sundstrom  suggest 
that  the  boundary  conditions  should  be  based  upon  their  ability  to  control  the  growth 
of  gravity  wave  disturbance  rather  than  on  their  accuracy  in  describing  the  quasi- 
geostrophic  part  of  the  motion.  In  this  connection  they  propose  that  only  those 
boundary  conditions  that  are  required  by  the  corresponding  linearized  system  of 
differential  equations  should  be  specified  a priori,  and  that  the  remaining  (computa- 
tional) boundary  conditions  should  be  obtained  from  the  interior  (and  boundary) 

values  of  the  limited  area  by  some  suitable  interpolation  (extrapolation)  procedures 

2 1 

This  is  essentially  the  approach  that  was  first  used  bv  Charney  et  al,  “ however,  as 

2° 

was  subsequently  shown  by  Platzman,  " they  used  an  inappropriate  (unstable)  extra- 
polation procedure  for  the  outflow  boundary  points.  To  avoid  the  instability  they 
abandoned  the  use  of  the  'proper'  boundary  conditions  and  resorted  to  simple  over- 
specification with  a smoothing  procedure.  For  the  primitive,  barotropic  shallow- 

23 

water  equations,  Charney  argued  on  the  basis  of  the  characteristic  quantities  for 
the  one -dimensional  problem  that  similar  boundary  conditions  are  required:  namely, 
the  normal  velocity  should  he  specified  at  all  boundary  points  and  the  potential  vorti- 
city  at  inflow  points. 

(Because  of  the  large  number  of  references  cited  above,  they  will  not  be  listed  here. 
See  Reference  Page  No.  67,  for  References  15  through  23.  ) 
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Sund strom  shows,  by  means  of  the  energy  method,  that  the  barotropic  vorti- 

citv  equation  in  a confined  region  is  properly  posed  with  boundary  conditions  of  the 
Charney-Fjortoft-von  Neumann21  type  (as  well  as  other  more  general  boundary  con- 
ditions) for  cases  with  open  boundaries  (that  is,  where  no  physical  boundaries  exist 
and  the  fluid  flows  in  and  out  of  the  boundary  region).  With  such  proper  boundary 
conditions,  the  solution  is  unique  and  convergent;  that  is,  a small  perturbation  on 
the  original  solution  will  remain  small  and  the  solution  is  continuously  dependent 
on  the  boundary  conditions  and  initial  data  as  long  as  certain  derivatives  of  the 
solution  remain  bounded.  Elvius  and  Sundstrom  propose  three  other  sets  of 


boundary  conditions  for  this  problem.  In  all  three  they  specify  the  tangential  velo- 
city at  inflow  boundary  points,  so  that  the  distinction  among  the  three  cases  con- 
cerns the  specification  of  an  additional  boundary  condition  for  all  boundary  points. 

This  second  boundary  condition  determines  the  manner  in  which  the  gravity  wave 
solutions  will  interact  with  the  boundary.  If  the  normal  velocity  is  prescribed  at 
the  boundary,  gravity  waves  are  essentially  trapped  and  reflected  back  into  the 
limited  area.  If  a combination  of  the  normal  velocity  and  the  geopotential,  corres- 
ponding to  the  ingoing  characteristic  for  the  gravity  waves,  is  specified  then  a 
gravity  wave  approaching  the  boundary  will  pass  out  of  the  region  without  reflection 
and  will  be  harmless.  If  on  the  other  hand  the  geopotential  itself  is  specified  for 
the  second  boundary  condition,  this,  like  the  first  case,  produces  reflection  of 
gravity  waves.  But  while  stable,  though  noisy,  difference  schemes  have  been  de- 
vised for  the  first  case,  no  stable  difference  scheme  is  known  for  the  third  case. 

The  number  of  additional  (computational)  boundary  conditions  (beyond  those 
specified  a priori)  depends  on  the  order  of  the  difference  equations.  The  additional 
conditions  should  not  be  specified,  but  should  be  determined  from  the  boundary  and 

interior  values  of  the  limited  area  bv  a method  that  is  both  accurate  and  stable. 

14 

In  the  final  analysis  Elvius  and  Sundstrom  assert  that  there  is  no  easy,  direct 
way  of  determining  the  stability  properties  of  two-  and  three-dimensional  differ- 
ence schemes  for  limited  area  problems.  The  energy  method  has  been  used  only 
for  the  case  with  boundary  conditions  that  are  open  for  external  gravity  waves.  For 

other  types  of  boundary  conditions  a method  of  stability  analysis  for  initial -boundary 

25  ,13 

value  problems  has  been  developed  by  Kreiss“  anil  applied  bv  (lustafsson  et  al, 
but  only  for  the  one -dimensional  case.  Consequently,  for  any  realistic  problem, 
we  are  forced  to  rely  on  computer  experimentation.  In  just  such  experiments, 

Elvius  and  Sundstrom1  have  obtained  reasonable  results  with  both  the  first  and 
second  types  of  boundary  conditions;  that  is,  where  the  tangential  velocity  is  speci- 
fied at  inflow  points  and  where  either  the  normal  velocity  or  a combination  of  normal 
velocity  and  geopotential  is  specified  at  outriow  points.  However,  in  those  experiments 

24.  Sundstrom,  A.  (1969)  Stability  theorems  for  the  barotropic  vorticitv  equation. 

Mon.  Wea.  Kev.  . '>7:340-345. 
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where  the  initial  data  contains  short  as  well  as  long  waves,  the  results  tended  to  be 
noisy.  The  primitive  barotropic  system  of  (shallow  water)  equations  is  a simple 
hyperbolic  system  which,  by  the  use  of  the  energy  method,  can  easily  be  shown  to 
be  well-posed  with  such  boundary  conditions.  However,  the  more  realistic  primitive 
system  of  equations,  which  describes  a baroclinic,  inviscid,  dry,  adiabatic  and 
hydrostatic  fluid,  does  not  form  a simple  hyperbolic  set  of  equations.  Although  the 
energy  transfer  equation  has  not  been  obtained  for  the  linearized  disturbance  equa- 
tions corresponding  to  this  system,  it  has  been  obtained  for  a simplified  system  for 
which  it  can  be  shown  that  the  pure  initial  value  problem  is  well  posed  in  that  a 
small  change  in  the  initial  conditions  will  remain  small  for  all  time  as  measured  by 
a RMS  norm.  ^ 

For  the  finite -difference  formulation  any  difference  between  the  prescribed 
boundary  conditions  and  the  correct  boundary  conditions  will  generate  gravity  waves 
and,  although  the  errors  may  remain  bounded,  the  solutions  will  be  incorrect.  If 
the  differences  between  the  prescribed  and  correct  boundary  conditions  are  large, 
nonlinear  effects  may  generate  real  disturbances  which  could  grow  rapidly.  A 
correct  set  of  boundary  conditions  should  specify  the  normal  velocity  and  tempera- 
ture in  such  a way  that  the  outward  propagating  components  are  not  hindered  in 
passing  out  of  the  region.  While  it  might  be  practicable  to  do  this  by  using  the  phase 
speed  of  a single  principal  component  in  a Sommerfeld  radiation  boundary  condition 

(as  suggested  by  Pearson ^®),  it  is  not  practicable  for  a multilevel  model  possessing 

27 

several  eigenfunctions.  To  overcome  this  difficulty  Orlanski”  suggests  a modifica- 
tion of  the  Sommerfeld  radiation  condition  in  which  rather  than  a single  constant 
phase  velocity,  a separate  phase  velocity  is  obtained  for  each  variable  from  the  grid 
point  values  in  the  vicinity  of  the  boundary.  This  approach,  which  involves  a cer- 
tain amount  of  inherent  smoothing  has  not  been  adequately  tested.  Instead  Sundstrom*' 
proposes  to  compromise  the  correctness  of  the  boundary  conditions  for  the  multilevel 
model  by  using  boundary  conditions  that  are  appropriate  for  the  shallow  water  equa- 
tion, that  is,  to  prescribe  the  tangential  velocity  at  inflow  points  and  compute  it  at 
outflow  points  by  a suitable  extrapolation  from  the  interior.  In  addition,  the  combina- 
tion of  the  normal  velocity  and  temperature  corresponding  to  the  first  inward  ex- 
ternal gravity  wave  is  prescribed  everywhere  and  the  combination  corresponding  to 
the  first  outward  external  gravity  wave  is  computed  by  extrapolation.  These  are 
approximately  the  barotropic  parts  of  the  normal  velocity  and  temperature.  Finally 
the  remaining  baroclinic  part  of  the  normal  velocity  and  temperature  is  specified  at 
all  boundary  points.  These  conditions  are  admittedly  an  overspecification  which 
permits  the  propagation  and  reflection  of  internal  gravity -wave  disturbances  into 


26.  Pearson,  R.  A.  (1974)  Consistent  boundary  conditions  for  numerical  models  of 

systems  that  admit  dispersive  waves,  J,  Atmos.  Sci.  , 31:1481-1489. 
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the  region.  Sundstrom18  satisfies  himself  with  the  hope  that  the  difficulties 
created  by  the  internal  gravity  waves  will  be  much  less  than  those  which  could  be 
generated  by  external  gravity  waves  and  that  their  effect  can  be  diminished  further 
by  adding  a diffusion  term  to  the  equation  near  the  boundary  zone. 

It  is  apparent  that,  even  if  it  were  possible  to  define  the  "correct"  boundary 
conditions  for  a realistic  baroclinic  limited-area  problem,  in  practice  the  boundary 
conditions  would  at  the  very  least  contain  truncation  errors  inasmuch  as  they 
would  be  determined  not  from  continuous  fields  but  from  a small  number  of  point 
values.  However,  in  reality  the  boundary  errors  will  be  more  severe  in  the  finite- 
difference  solution  since  the  proper  boundary  conditions  cannot  be  defined  for  realis- 
tic problems  and  at  test  one  can  apply  boundary  conditions  which  are  appropriate 
for  much  simpler  problems.  Nevertheless,  even  for  the  simpler  shallow- water 
problem,  obtaining  the  proper  boundary  conditions  is  likely  to  involve  computa- 
tional and  programming  complexities.  In  this  connection,  Williamson  and  Brown- 
28 

ing  use  a relatively  simple  set  of  lateral  boundary  conditions  for  the  NCAR 

2Q 

limited-area  model.  Following  the  approach  of  Shapiro  and  O' Brien  ' who  worked 

with  the  nondivergent  barotropic  vortieity  equation,  Williamson  and  Browning28 

specify  u,  v,  and  t on  points  determined  as  inflow  points  from  the  global  model 

solution.  At  outflow  points  they  determine  these  parameters  by  interpolation  from 

the  fine-mesh  solution,  bv  determining  the  interior  point  which  would  have  been  the 

starting  point  of  a parcel  following  a trajectory  determined  from  the  coarse-mesh 

flow  field.  As  would  be  expected  their  results  show  appreciable  two-grid-interval 

noise,  but  the  noise  does  not  appear  to  interfere  with  the  definition  of  the  large- 

scale  fields.  Nevertheless,  they  find  it  expedient  to  apply  additional  smoothing  in 

a band  of  several  grid  intervals  around  the  boundary  of  the  limited  area. 

30  31 

Davies  ’ makes  use  of  the  energy  method  to  establish  criteria  of  uniqueness, 
which  he  then  uses  to  determine  boundary  conditions  for  the  barotropic  shallow  water 

equations  as  well  as  for  the  baroclinic  primitive  equations.  However,  these  condi- 

3 2 

tions  have  been  criticized  for  overspecification  by  DeRivas  as  well  as  bv 
..15  33 

Sundstrom.  Recently,  Davies'  has  proposed  a different  procedure.  Starting 

28.  Williamson,  D.  L. , and  Browning.  G.  L.  (1974)  Formulations  of  the  lateral 

boundary  conditions  for  the  NCAR  limited-area  model,  .1.  Appl.  Meteor.,  I 

13:8-16.  
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29.  Shapiro,  M.A.,  and  O'Brien,  (1970)  Boundary  conditions  for  fine-mesh 

limited-area  forecasts,  .1,  Appl.  Meteor.  , j):345-349. 

30.  Davies,  II.  C.  (1973a)  On  the  lateral  boundary  conditions  for  the  primitive 

equations,  .1.  Atmos.  Sci.  , 30:147-150. 
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from  the  position  that  in  all  practical  applications  (in  the  one-way  interaction 
limited  area  problem)  the  available  boundary  data  will  contain  inherent  errors,  lie 
proposes  an  expedient  method  involving  the  relaxation  of  the  interior  flow  in  the 
vicinity  of  the  boundary  to  the  external  fully  prescribed  flow.  The  method  reduces 
the  sensitivity  of  the  governing  equations  to  overspecifications  of  the  boundary  data 
ind  achieves  a reduction  in  the  noise  generation  at  the  lateral  boundaries  by  means 
of  an  artificial  modification  of  the  equations.  The  modification  consists  of  a 
Newtonian-type  adjustment  term  with  a relaxation  coefficient  whose  spatial  depen- 
[ dence  is  such  as  to  force  the  interior  field  values  to  approach  the  specified  boundarv 

values  in  a smooth  fashion.  The  method  has  the  virtue  of  simplicity  and  seems  to 
be  effective  for  the  single  test  cases  shown;  however,  it  is  not  obvious  to  what 
physical  problem  the  adulterated  equations  apply.  While  it  is  apparent  that  the 
method  induces  a certain  amount  of  smoothing  between  the  specified  boundarv  values 
and  the  interior  fine-mesh  values  tiie  total  smoothing  function  does  not  have  a single 
amplitude  response  function  (because  of  the  spatial  variation  in  the  relaxation 
coefficient)  and  therefore  it  is  not  possible  to  assess  the  extent  of  the  smoothing 
involved. 

That  some  form  of  smoothing  might  be  necessary  can  be  inferred  from  the  re- 
cent comments  of  Dennett.  ' * lie  asserts  that  suitable  open  boundary  conditions  at 

I the  outflow  points  of  a limited  area  (that  is,  boundary  conditions  which  avoid  trap- 

ping computational  noise  in  the  limited  forecast  area)  are  possible  only  for  small 
amplitude,  short  wavelength,  barotropic  waves;  that  noise  problems  are  unavoid- 
able in  planetary  or  baroclinic  wave  forecasts.  " 

On  the  assumption  that  for  realistic,  three-dimensional  problems  a certain 
amount  of  boundary  error  is  inevitable  (regardless  of  the  choice  of  boundary  condi- 
tions and  the  method  used  to  obtain  them)  a high-order  interpolation  procedure  was 
proposed'^1  for  supplying  time-dependent  boundarv  values  in  a fine-mesh  limited- 
area  model.  The  interpolation  was  performed  on  a ’global1  coarse- mesh  solution 
which  was  run  concurrently  with,  but  independent  of,  the  fine-mesh  solution,  l’he 
aim  of  this  experiment  was  to  determine  whether  the  use  of  a high-order  inter- 
polation (eight  points  in  one  dimension',  which  was  designed  to  represent  faithfully 
the  amplitudes  of  all  waves  exrept  the  very  short  waves,  ,f  would  produce  a smooth 
transition  between  the  specified  boundary  values  and  the  calculated,  interior  fine- 

mesh  values.  Numerical  simulation  experiments  were  carried  out  in  a nested 
TT  Bennett,  A.p.  (107a)  Open  Boundarv  Conditions  for  the  Limited  Area  I orc- 
c asting  Problem , Report  Vo.  0.  The  A RB  Programme  on  Numerical 
Pxpetu mentation,  A.  Robert,  Pd.,  pp  114-11.1 

•if).  Shapiro.  H.  Hn7:i>  \ I liu li  -Or  ler  Interpolation  Procedure  lor  I se  in  rinc-Mi  s' 
lamited-Area  Modi  Is , Phvs.  Sei.  11  e s . 1 ’up.  , No.  I"' fid,  \ It-  R I -TR  -t  1 -dl-t  •' . 
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region  with  a grid  spacing  one-fifth  the  size  of  the  coarse  grid  length  within  the 
framework  of  an  eight-level,  drv,  inviscid,  primitive-equation  model.  Boundary 
information  was  supplied  from  the  coarse-mesh  solution  by  interpolation  for  all 
prognostic  variables  for  all  of  the  fine-mesh  boundary  points,  regardless  of 
whether  the  flow  was  in  or  out  of  the  limited  area.  This  is  obviously  an  over- 
specification of  the  boundary  conditions,  but,  as  shown  by  Chen  and  Mivakoda,  1 
such  overspecification  has  the  advantage  of  simplicity  and  with  smoothing  gives 
results  comparable  to  those  obtained  with  presumably  well-posed  boundary  condi- 
tions. 

35 

Tn  the  above-mentioned  study  a comparison  was  made  between  two  methods 
of  interpolation  from  the  coarse-mesh  solution,  a linear  and  a high-order  interpola- 
tion. The  high-order  interpolation  procedure  is  strongly  damping  for  two-  and 

three-grid  interval  waves,  but  barely  has  any  effect  on  waves  longer  than  about  six 

3 6 

grid  lengths  even  after  many  applications  of  the  operator.  The  linear  interpola- 
tion is  at  least  as  strongly  damping  for  short  waves  as  the  high-order  interpolation, 
but  more  importantly,  the  linear  interpolation  is  strongly  damping  for  even  moder- 
ately long  waves. 

The  solutions  obtained  with  the  use  of  high-order  interpolation,  although  by  no 
means  noise-free,  were  noticeably  smoother  than  those  obtained  with  linear  inter- 
polation. The  latter  solutions  were  seriously  contaminated  with  noise.  Since  there 
was  no  viscous  dissipation  in  the  model,  a weak  smoothing  operation  was  applied 
to  the  coarse-mesh  solution  every  ten  time  steps  (every  200  min).  This  same 
smoothing  was  also  applied  to  each  fine-mesh  solution  every  200  min  (every  50  time 
steps).  In  view  of  the  weakness  of  the  smoothing  operation,  the  relative  smoothness 
of  the  solution  with  high-order  interpolation  was  striking.  It  was  felt  that  the  re- 
maining roughness  of  these  solutions  could  be  controlled  by  a combination  of  more 
frequent  smoothing  and  by  using  an  interpolation  procedure  which  not  only  is  efficient 
in  restoring  wave  amplitudes  damped  by  linear  interpolation,  but  also  minimizes 
phase  error.  The  present  study  describes  an  attempt  to  remove  the  remaining 
roughness  by  such  means,  using  the  same  basic  model  framework  and  maintaining 
the  same  simple  overspecified  boundary  conditions. 

A modification  of  the  high-order  interpolation  procedure  used  in  the  earlier 
study  was  made  which,  while  still  maintaining  the  amplitude-restoring  properties, 
minimizes  the  phase  error  arising  from  interpolation.  A new  series  of  fine-mesh 
limited- area  solutions  is  carried  out  in  which  comparisons  are  made  among  three 
methods  of  interpolating  the  coarse-mesh  solution  onto  the  limited-area  boundaries: 
linear  interpolation,  high-order  amplitude  restoring  interpolation,  and  high-order 
phase  and  amplitude  restoring  interpolation.  In  addition,  solutions  are  obtained 
with  differing  frequencies  of  smoothing  in  order  to  test  the  relative  importance  of 

37.  Chen,  .1.11.,  and  Mivakoda,  K.  (1974)  A nested  grid  computation  for  the  baro- 
tropic  free  surface  atmosphere,  Mon.  Wen.  Rev.,  102.181-190. 
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the  smoothing  and  interpolation  operations.  The  aim  of  these  experiments  is  the 
development  of  a simple,  practical  procedure  for  solving  a multilevel  primitive 
equation  model  over  a limited  area  wit!)  fine  mesh  resolution  without  the  need  for 
solving  the  global  model  with  a fine-mesh  grid.  We  therefore  propose  to  specify 
all  prognostic  variables  at  all  boundary  points.  We  justify  this  overspecification 
of  the  boundary  conditions  on  the  basis  that: 

(11  The  proper  boundary  conditions  can  be  shown  only  for  simplified  linear 
systems. 

(2)  liven  if  wo  knew  the  proper  boundary  conditions  for  realistic,  three- 
dimensional  nonlinear  problems,  we  would  have  to  specify  some  of  the  boundary 
conditions  in  the  finite-difference  solution  from  the  outside  region,  where  at  best 
the  values  are  known  less  accurately.  C onsequently,  boundary  errors  are  in- 
evitable in  practice. 

(Ml  Since  we  do  not  know  the  proper  boundary  conditions  and  cannot  eliminate 
boundary  errors  completely,  some  kind  of  artificial  smoothing  is  necessary. 

(41  In  view  of  the  availability  of  simple  and  highly  effective  smoothing  and 
filtering  procedures  which  can  be  designed  to  be  as  scale-dependent  as  desired,  ,l<' 
and  in  view  of  the  results  of  Chen  and  l\1 iyakoda'*"  with  overspecified  boundary  con- 
ditions, it  is  felt  worthwhile  to  avoid  the  complexities  of  attempting  to  approximate 
proper  boundary  conditions  and  adapt  the  simple  expedient  of  overspecification  and 
filtering. 


2.  IVIT'.RI’OI.  \TI()> 

The  process  of  interpolation  is  widely  vised  to  obtain  estimates  of  parameters 
at  certain  points  in  time  or  space  from  known  or  measured  values  which  are  avail- 
able only  at  other  points.  It  is  tacitly  assumed,  when  using  interpolation,  that  the 
parameter  has  a value  at  the  point  in  question,  it  is  also  assumed  that  the  value  of 
the  parameter  at  any  desired  point  may  he  estimated  from  some  function  of  the 
known  or  given  values  of  the  parameter.  In  effect,  it  is  assumed  that  the  sample 
of  the  known  set  of  values  is  representative  of  the  entire  field.  If  we  make  explicit 
what  is  usually  implicit  in  anv  event,  namelv  that  the  parameter  is  continuous  over 
tiic  domain,  we  may  speak  of  its  Courier  spectrum.  In  particular,  we  consider  the 
discrete,  finite  Courier  representation  of  the  field  which  is  obtained  from  the  set  of 
values  at  known  points  to  contain  all  of  the  available  information  concerning  the  field 
and  that  the  best  estimate  of  the  value  of  the  parameter  if  tile  point  in  question  is 
that  which  is  derived  from  this  Courier  representation  of  the  known  field.  1‘sing 

38.  Shapiro,  H.  ( 10701  Smoothing,  filtering,  .nd  boundary  effects,  I lev,  (leophvs. 
and  Space  I’hvs.  , 8:350-387. 
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this  measure  of  best  estimate,  the  author  proposed  an  interpolation  operator 
which  compensates  for  the  Fourier  amplitude  truncation  produced  by  two-point 
linear  interpolation.  The  aim  of  this  report  is  to  generalize  the  operator  and  to 
propose  a specialized  procedure  which  minimizes  phase  as  well  as  amplitude  trun- 
cation. 


2.1  Amplitude  Restoration 

39 

Following  the  procedure  adopted  earlier,  ‘ " let  fix)  represent  some  function  of 

time  or  space  in  one  dimension,  integr.ible  in  the  finite  domain  -D  < x < D,  and  let 

fix  ) = f.  for  discrete  values  of  x such  that  x.  = iAx,  where  i is  an  integer  and 
11  i 

Ax  0 is  such  that  2D/Ax  is  an  integer.  Two-point  linear  interpolation  at  the  point 
i = rAx,  where  0 £ r s 1,  is  defined  by 


-(0)  . 
f.  = rf.  , 
i+  r i+l 


<1  - r>  f. 


(1) 


where  the  interval  (-1),  D)  has  been  scaled  so  that  Ax  = 1. 

Then  f.  can  be  expressed  in  terms  of  a sum  of  Fourier  components  of  the 

general  form  A sin  (nx.  + t?  ),  where  A is  the  amplitude  of  the  wave  component 
* n i n n 

with  wave  number  n in  = 2 i/>,  where  A is  the  wavelength  of  the  component',  and 
is  the  phase  of  the  component;  f the  value  of  the  function  at  the  point  i+r  can 
be  similarlv  represented  bv  a sum  of  components  of  the  form  A sin  I nix.  + rAx)  + o 

(0)  n L 1 "J 

and  the  corresponding  interpolated  value,  f ^ by  components  of  the  form 


(0) 


A<°)  _ 

n 


in  nix.  + rAx)  + O - d , where 
L i n nj 

i / 1 , . 2 nAx  I 1 “ , 

4 r ( 1 - r)  sin  — « — A , 

2 J n 


(2) 


and  where 


-1  i (1-r)  sin  nrAx  - r sin  [nil  - r)  Ax]  I _ 

n ~ tan  | (1-r)  cos  nrAx  + r cos  ] n ( 1 - r)  Ax]  j ’ 

A^0'  A ^ expresses  the  amplitude  damping,  and  d the  ptiase  stiift,  produced 

by  two-point  linear  interpolation,  it  is  apparent  that  two-point  linear  interpolation 

can  introduce  considerable  damping  of  amplitude  and  shifting  of  phase,  especially 

for  the  higher  wave  number  components.  For  wave  numbers  greater  than  zero 

(omitting  the  cases  where  r - Oor  I for  which  f*^'  - f. , Kd  0 only  if  r 0.  f>.  it 

& i+r  i+r 

is  also  apparent  from  (3)  that  dir)  = - d (1  - r). 

The  value  of  r (or  1 - r'  for  which  the  magnitude  of  d is  a maximum  is  a function 
of  n and  can  easily  be  determined  from  (3).  For  waves  equal  to  or  longer  than  3-grid 
intervals,  this  value  of  r is  near  0.  2!i.  It  varies  from  0.  213 f> f>  for  3-grid  interval 
waves  to  0.  2113  for  infinitely  long  waves.  For  the  2-grid  lntrrv  d wave,  Fq.  (3) 
reduces  to  d a nr,  except  for  r 0,  1/2,  or  I,  for  which  d 0.  Thus,  for  the 


■ 
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2-grid  interval  wave,  there  is  a discontinuity  in  d at  r = 0.  5.  For  all  longer  waves, 
d is  continuous. 

It  has  been  shown***  that  the  inverse  of  for  all  values  of  t and  a (except  for 
9 n 
t sin"  a = 1 for  which  there  is  no  inverse)  is  given  by  the  infinite  series 


(0)1  -1  , t . 2 1.  3 .2  . 4 1.  3. 5 ,3  6 

p ^ — 1 + — sin  a + ^ t sin  a + ^ — jr  t sin  a 


+ . . . + 


1.  3.  5.  ...  <2h  -1)  Ji  . 2h 


2.t;-(5.  — m) 


t sin  a + . . . , 


(4) 


where  t = 4r  (1  - r)  and  a = nAx/2.  Then  t sin"  a = 1 occurs  onlv  for  the  2Ax  wave- 
length with  r = 0.  5.  This  wavelength,  the  shortest  Fourier  component  that  can  be 
resolved  with  the  given  data,  is  completely  eliminated  by  two-point  linear  interpola- 
tion when  r = 0.  5. 

The  series  (4)  is  a monotonicallv  increasing  function  of  h.  If  we  operate  on 

f|^|_  with  an  operator  representing  a truncation  of  the  series  (4),  we  can  (with  the 

exception  noted  above)  restore  the  amplitudes  of  the  various  Fourier  components 

2 

which  had  been  damped  by  two-point  linear  interpolation.  We  note  that  sin"  a is  the 

2 

Fourier  representation  of  the  operator  (-6  / 4 ) where  6(  ).  = ( 'i-l/9‘ 

That  is,  operating  on  any  function  with  (-62/4)  is  equivalent  to  multiplying  the  ampli- 

2 

tude  of  the  component  with  wave  number  n bv  sin  a. 


6 ( ).  = ( 


).  . -2  < 
i+  1 


).  + ( ).  , . 
i i - 1 


(5) 


Thus,  the  series  (4)  can  be  represented  as  the  operator: 


(o7 

_ i 

( -t  .2  \ __  1.3 

/ "»  A2 

9 

\ ” 

1,  3.  5 . . 

. (2h- 1)  / 

Lp-j 

= 14  1 /2  | 

T6  ) 4 77T 

( T6 

/ + •• 

• + Trmrrr 

. " /2FiT  \ 

(4«2) 


(<;> 


It  is  apparent  that  the  operator  defined  bv  (fi)  operating  on  fj^J.  is  linear  (additive),  in 
that  each  term  of  the  series  consists  of  a centered  symmetrical  operation  involving 
two  grid  points  additional  to  the  next  lower  ordered  term.  Thus,  the  interpolation 
operator  of  2 (h  + 1)  points  in  one  dimension  consists  of  the  operator  defined  by  the 
series  (G)  truncated  at  the  hth  term,  operating  upon  fj^J,.  This  operator  differs  from 
the  operator  discussed  in  Shapiro**’  which  is  multiplicative  and  less  efficient.  It  is 
apparent  that  this  operator,  which  we  can  write  as 


,( h) 


1+ 


Z\  1.3/-t,2\  1.3,5.  ...  (2h-l)  Hr2\ 

1/2  ( T 6 ) 4 77T  I T 6 ) 4 ?.4.  R.  ...  (?!,)  It6/ 


,<0) 
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h is  an  amplitude  response  which,  in  the  limit  as  h approaches  infinity,  approaches 
2-1/2  21/2 

(1  - t sin"  a)  ( 1 - t sin  a)  = 1.  Furthermore,  we  shall  show,  for  any  value 
of  h,  the  amplitude  response  is  a maximum. 

For  am  level  of  restoration  h,  and  anv  function  F. 


„2h 


tion  of  (2h  + DF^,  k = ,j,  j ± 1, 


C2hF..(-ll1,'2h 


h- 1 
T (-1) 
k=0 


, j ± h,  given  by 
k / 2 


.1 


F.  is  a linear  combina- 

.1 


(?) 


* j+(h-k)  4 * j-(h-k) 


(7) 


where  (^1  are  the  binomial  coefficients.  Thus  (6)  truncated  at  the  level  h,  may  be 
written  as  an  operator  on  F. 


R(h)  (F .)  = F.  + 1 

.1  .1 


2F . + F . . + F.  , 
.1  j+1  1-1 


2 r 

1 6F.  - 4 

f F- 

+ F , 

+ F. 

+ F. 

) [_  j 

V .1+1 

j-1 

3+2 

j -21 

-J 

. 3.  5 
TTH? 

C<f[ 

-20F. + 15  / 

1 

F.  , + F.  , 
J+l  I"  1 j 

1 "6  (’>2*  Fi-2)  * 

v.]  * 

1.  3.  5. 

...  (2h-l) 

uf 

j+3  + 

• • * 4 2.4.6. 

...  (2h) 

h-1  2h 

v <-l>k  k | ( F. 
k=0  ' 


i+(h-k)  * 1 j-(h-k) 


(8) 


If  F.  represents  two-point  linear  interpolation,  f!^\  then  K^1*  (F.)  = ff*1'  is  an 
) ^ r i+r  j l+r 

interpolation  operator  which,  at  the  level  of  restoration,  h,  maximizes  the  restora- 
tion of  amplitude  damped  bv  f((".  That  is,  the  interpolation  operator  defined  by  (8' 
(01  i+r 

operating  on  f r is  ideal  in  the  sense  that  for  any  linear  operator  involving  2(h  + 11 
points  it  maximizes  the  restoration  of  amplitude  damped  by  two-point  linear  inter- 
polation. 

The  property  of  maximum  restoration  can  be  demonstrated  by  induction.  Using 
Eqs.  (2)  and  (4),  we  write  the  amplitude  response  function  for  the  hth  level  of 
restoration  as: 


(hi  (oi  r 

, t .2  3 2 .4 

Pn  = Pn 

l+TjSin  a + -j  t sin  a 

. . + 

1.3.5.  , 

. . . ( 2h  - 1 1 .h  . 2h  1 

4 2.4.6. 

1 sm  a 

J 

Then  we  write  amplitude  response  function  for  the  (h  + 11th  level  of  restoration  as 
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1 


P(h+1)  = p<0) 


, t . 2 1.  3. 5.  ...  (2  h-1)  ,h  . 2h 

1 + -g-  sin  a 4 . . 4 — JTTT! — t sln  a + 


,.h4l  . 2(h 

4 At  sin 


1) 

a 


We  solve  for  A under  the  condition  that  dp  ' */dh  = 0 and  p 1 * is  a maximum 

' n f n 

and  find  that 


A = 


1.  3.  5. 


(2  h 4 1) 


7 4757  . 777  + 2T 


which  is  the  coefficient  of  t*1+  * sin"^+ " a in  the  series  (4). 
Inasmuch  as  the  general  term  in  (8) 


1.  3.  5.  . . . (2h  - 1) 

77475."—“  72hl 

(2h)  ! 


(2h)  I 


[2h.  (2h  - 2). 
2h  ' 


2^h: 


T ( H ) • 
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with  the  aid  of  (7),  (8)  can  be  expressed  in  compact  form 


R(h)  (F.)  = T 

' h=  0 


(t/4) 

TFT 


2 

)4(h-k)  + 1 j-(h-k)J 


2h 

h 


F.  4 

.1 


<-ir<t/4> 


7717 


>(?) 


h-l  , 
v (-1)" 
k-0 


(IN 


(9) 


Similarly,  the  amplitude  response  function  corresponding  to  the  amplitude- restoring 
interpolation  operator  (9)  is 


P<h>  = p<0)  y 
h 


n ‘ n 
JO)  .. 


(?) 


..,,.h  . 2h 

(t/4)  sin  a 


where  is  the  amplitude  response  function  corresponding  to  the  operator  F.  in 
(9).  In  all  applications  in  this  paper  F^  represents  two-point  linear  interpolation 
at  the  point  j = i + r.  Thus,  R**''  (F.)  = f|^*|,  in  any  application. 

Figure  I demonstrates  the  etticiencv  of  (9>  as  an  amplitude  restoring  interpola- 
tion operator  for  various  values  of  h for  the  case  where  r - 0.  5,  for  which  the  ampli- 
tude error  is  a maximum  for  any  given  h.  The  values  of  h range  from  0,  which  cor- 
responds to  two-point  linear  interpolation,  to  8,  which  corresponds  to  an  18-point 
interpolation  operator.  The  amplitude  response  of  each  operator,  represented  as 
P in  Figure  1,  is  shown  as  a function  of  wave  number  determined  from  7 2 equally 
spaced  grid  points.  Wave  number  36,  therefore,  corresponds  to  a wavelength  of 
2Ax,  twice  the  grid  distance.  Since  the  two-grid  interval  wave  is  eliminated  by 
two-point  linear  interpolation  when  r = 0.  f>,  the  amplitude  of  this  wave  remains  zero 
for  all  h.  The  ordinate  scale  in  Figure  1 is  linear  from  0.  0 to  0.  9 and  logarithmic 
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above  0.  9.  Since  fractional  wave  numbers  have  no  meaning,  the  curves  connecting 
the  points  for  each  operator  have  been  drawn  to  assist  in  the  interpretation  of  the 
figure.  The  curves  are  not  joined  for  points  on  opposite  sides  of  the  0.  9 dividing 
line  to  emphasize  the  change  of  scale. 


Figure  1.  The  Amplitude  Response  of  the  Amplitude  Restoring  Interpolation 
Operator,  Kq.  (91,  as  a Function  of  Wave  Number  in  a 72-I'oint  Periodic 
domain.  The  operator  is  represented  as  P^  where  h indicates  the  order  of 

restoration.  The  ordinate  scale  is  linear  from  0.  0 to  0.  9 and  logarithmic  i 

for  higher  values 

It  is  apparent  from  Figure  1 that  there  is  almost  complete  restoration  of  the 
amplitudes  of  the  longer  waves  even  for  relatively  low  values  of  h.  For  values  of 
h > 3,  there  is  very  little  damping  of  waves  longer  than  3-grid  intervals  (wave  number 
241.  For  example,  the  amplitude  of  the  3-grid-interval  wave  for  li  4 is  0.  902,  as 
compared  with  its  value  for  h = 0 of  0.  900,  and  for  the  4-grid-interval  wave  (wave 
number  181,  the  comparable  values  are  0.990  and  0.707.  Since  the  damping  is  a 
maximum  for  r = 0.  5,  the  restoration  of  amplitude  for  any  h is  somewhat  more 
complete  for  other  values  of  r than  that  shown  in  Figure  1.  For  any  given  r,  the 
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restoration  of  amplitude  is  a function  only  of  h and  wave  number;  consequently,  the 
most  suitable  level  of  restoration  can  be  selected  for  a particular  application.  How- 
ever the  ability  of  an  interpolation  operator  to  approximate  a given  function  depends 
upon  phase  error  as  well  as  amplitude  error.  Therefore  we  shall  consider  the 
merits  of  an  interpolation  operator  which  minimizes  both  amplitude  and  ph  :se  error. 

2.2  lliasc  Kti'sloiailion 

As  we  have  seen,  two-point  linear  interpolation  not  onlv  T mps  the  amplitude 
of  e ich  l ourier  component  but  also,  in  general,  shift-  the  phase  of  each  component, 
f igure  2 shows  the  phase  difference,  d,  1 q.  <3>]  in  radians  is  a function  of  r,  is 
well  as  -d  as  a function  of  I 1 - r'  in  parentheses,  for  waves  longer  than  3-grid 
intervals.  1'he  onlv  exception  to  phase  shifting  occurs  when  the  interpolation  point 
is  midwav  between  ttie  two  data  points;  that  is.  when  r - 0.  f>.  This  feature  of  two- 
point  linear  interpolation  for  r 0.  a permits  the  definition  of  a special  interpolation 
procedure  for  the  case  where  the  r itio  of  the  original  grid  length  to  the  final  inter- 
polation grid  length  is  ..  power  of  2. 


\\ 


v.  1 4*1 


figure  2.  1’he  Phase  Difference,  d or  (-d),  in  Hadians  for  Two-Point 
I .inear  Interpolation  as  a function  of  Wavelength  in  C'.rid  Intervals  for 
Various  Values  of  r or  (1-rl 


If  we  apply  two-point  linear  interpolation  | fq.  ( PI  with  r 0.  !i  throughout  the 
domain,  we  obtain  values  of  the  parameter  at  a set  of  discrete  points  consisting  of 
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original  data  points  and  interpolated  data  points  situated  midway  between  the  original 
points.  Since  no  phase  shift  lias  been  introduced  by  this  interpolation  procedure,  we 
may  correct  for  the  amplitude  truncation  by  applying  the  high  ordered  amplitude 
restoring  operator  (hi  at  some  suitable  level  of  truncation,  say  h = q.  After  apply- 
ing (f*l,  the  set  of  original  and  interpolated  data  has  virtually  the  same  spectral 
properties  as  the  original  data  except  for  the  lisappearance  of  the  2 -grid-interval 
wavelength.  Since  the  true  2Ax  wave  is  likely  to  be  poorly  represented  in  the  orig- 
inal lata  and  troublesome  in  computations  because  it  is  often  strongly  contaminated 
bv  noise,  it  is  frequently  advantageous  to  eliminate  it  from  the  original  data.  There- 
fore, the  suppression  of  the  2Ax  wave  is  a desirable  and  advantageous  feature  of  the 
interpolation  operator.  1'he  set  of  combined  original  data  and  interpolated  data  w ith 
restored  amplitude  may  be  considered  as  a new  set  of  homogeneous  data  with  values 
given  at  grid  intervals  Ax  2.  The  above  procedure  can  now  be  repeated,  to  obtain 
in  turn,  other  sets  of  homogeneous  data  with  values  at  grid  intervals  Ax/4,  Ax  '8,  . . . 
Ax  '2m.  This  process  can  be  expressed  formally  as  follows; 

Up  define  a variable  position  point  a and  an  interim  grid  interval  Ax  such  that 

tv  - i • p Ax  21”  * anil  Ax  Ax/!l"  where  m 0,  1,2 M and  where  p 1, 

M * 1 ) 

3,  5,  ...  , (2  -1).  It  is  apparent  that  O'  represents  a special  class  of  position 

points  i i-  r*  between  i and  i • 1 where  r';  assumes  only  the  values  of  1 2 (for  n - 0); 
1/4  and  3/4  (for  m = 1);  1/8,  3/8,  5/8,  and  7/8  (for  m - 2);  and  so  forth.  It  is  also 
apparent  that  this  phase  and  amplitude  restoring  interpolation  operator  n av  be  ex- 
pressed in  terms  of  the  amplitude  restoring  operator  (!)),  as. 


f (h.q)  H(h)  j,  (q))  R(h)  , / 2 (f  ' . f <q*  ) . (101 

f(/’-q<  is  thus  seen  to  consist  of  an  amplitude  restoring  operation  at  level  h on  t he  opera  - 

tori’,  which  in  turn  consists  of  a two-point  linear  interpolation  on  the  fiel.i  iff  1,1  at 

<«>  Ax'*  Ax  1 

the  point  (w  hich  is  midway  between  the  interim  grid  points  <v  ‘ — n — and  o — . The 

f.tq*  consist  of  the  original  grid  point  values  of  the  function  < f . , f . _ . , ...  1 as  well  as 
interpolated  values  obtained  from  previous  applications  of  (10)  with  amplitude  re- 
storation at  the  level  q.  To  the  extent  that  the  f.'q'  are  homogeneous,  there  will  be 

no  phase  error  in  the  f Since  the  f/q*  consist  of  both  original  and  inter- 

« .1 

polated  grid  point  values  of  f(x),  they  cannot  be  completely  homogeneous,  but  the 
degree  of  homogeneity  depends  on  the  level  of  restoration  (q)  and  is  therefore  to  a 
large  extent  controllable.  After  final  application  of  (101,  with  m - \t,  original  grid 
point  values  of  the  function  f(x),  or  interpolated  values  with  corrected  phase  and 
amplitude,  are  available  at  all  points  i + s/2^*  * where  s - 1,  2,  3,  ...  , (2^'*  '-11. 

The  growth  of  the  size  of  the  stencil  of  operator  (101  depends  not  only  on  h and 
q but  also  on  m and  p.  Thus,  when  m - 0,  (101  reduces  to  the  amplitude  restoring 
operator  (0),  and  the  number  of  data  points  involved  in  the  stencil  represented  by 


( 10)  is  12  (h  * 1);  when  tn  1 , and  q h the  number  is  1?  3 h;  but  w lien  m is  equal  to 

or  [greater  than  2,  the  number  of  terms  also  depends  upon  p which  letermines  the 
position  of  rv  between  the  points  i and  i • 1.  For  example,  as  ran  he  seen  fron 
Table  1,  when  m 2 and  h 1,  there  are  5 terms  in  the  stencil  for  r - 1/fl,  but 
6 terms  for  r 3/8.  On  the  other  hand,  for  the  same  value  of  m,  but  with  h 3, 
there  are  13  terms  for  the  r - 1 ' 8 position  and  only  12  terms  for  the  r - 3/8  posi- 
tion. 


Table  I.  Coefficients  of  the  Interpolation  Stencils  for  \mplitude  and  Phase 
llestoration  in  Terms  of  the  Original  Grid  Point  Values  for  Various  Values 
of  m,  r,  and  h (where  q is  taken  equal  to  h) 


3 : m = 0 ; r 1/2 


h 0 

1 Hi  * 

h 2 

h = 3 

r - 

h 4 

Zi-4Ax 

. 00053406 

Zi-3Ax 

-.00244140 

-.00617881 

Zi-2A.x 

. 0117  1875 

. 02392578 

. 03460683 

Zi-Ax 

-.  0625 

-.  087  65625 

-.  11862881 

-.  13458252 

z. 

1 

. 5 

. 5625 

. 58583750 

. 58814453 

. 60562134 

ZiAx 

. 5 

. 5625 

. 58583750 

.58814453 

. 60562134 

Zi*2Ax 

0625 

08765625 

-11862881 

. 13458252 

Zi<-3Ax 

.0117  1875 

. 02382578 

. 03460683 

Zi*-4  Ax 

-.  00244140 

-.  006  1788  1 

z.  - 

U 5 x 

. 00053406 

B:  m-  1;  r = 1/4 


. 00380625 
. 07031250 


h = 2 

h = 3 

h = 4 

-.  00000330 

. 00000586 

-.  00003368 

. 00023365 

. 00084127 

-.00114441 

00403047 

-.  007  54884 

. 01640320 

. 027  10342 

. 0356480 1 

-.  10258484 

-.  12007284 

-.  13095344 

Table  1.  Coefficients  of  the  Interpolation  Stencils  for  Amplitude  and  Phase 
Restoration  in  Terms  of  the  Original  Grid  Point  Values  for  Various  Values 
of  m,  r,  and  h (where  q is  taken  equal  to  hi  (Conti 


B:  m = 1; 

r = 1 /4  (Conti 

h = 0 

h = 1 

h = 2 

h = 3 

h = 4 

z. 

1 

. 75 

. 84375000 

. 87089539 

. 88179588 

. 88708519 

z.  1 

1+ 1 

. 25 

. 25781250 

. 26206970 

. 26671171 

. 27098010 

03515625 

05149841 

06254196 

07102302 

Zi+3 

. 00572205 

. 01174092 

. 01745304 

Zi+4 

. 00013733 

. 00088787 

00268873 

Zi+5 

00005841 

. 00012805 

Zi+6 

. 00001518 

Z.  - 

1-f  < 

. 00000029 

C:  m = 2; 

r = 1/8 

h = 0 

h = 1 

h = 2 

h = 4 

Zi-7 

00000002 

Zi-6 

00000001 

00000382 

Zi-5 

. 00000998 

00002444 

Zi-4 

00001341 

. 00019600 

. 000697  24 

Zi-3 

00109041 

00303826 

00520480 

z,-2 

. 00439453 

. 01241863 

. 01854892 

. 02324251 

zi-l 

05175781 

06990135 

07  875438 

08409304 

zi 

. 875 

. 94921875 

. 96389651 

. 96858347 

. 97055880 

zi+l 

. 125 

. 1 1425781 

. 11541545 

. 1 1828768 

. 12100399 

Zi+2 

01611328 

02330482 

02866947 

03299573 

Zi+  3 

. 00251234 

. 005234  94 

. 00792706 

Zi+4 

. 00006706 

00036325 

00115760 

Zi+5 

00003006 

. 00004125 

Zi+6 

. 00000014 

. 00000854 

Zi+7 

. 00000007 

Zi+8 

+. 00000000 

I 
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Table  1.  Coefficients  of  the  Interpolation  Stencils  for  Amplitude  and  Phase 
Hestoration  in  Terms  of  the  Original  Grid  Point  Values  for  Various  Values 
of  m,  r,  and  h (where  q is  taken  equal  to  h)  (C'ont) 


As  h increases,  the  number  of  terms  in  the  stencils  increases  regardless  of 
the  value  of  m or  p,  but  for  h-  3 the  additional  terms  are  small  in  magnitude.  Thus, 
it  appears  that  as  far  as  phase  and  amplitude  errors  are  concerned,  for  ill  but  the 
shortest  wavelengths,  there  is  little  advantage  in  using  orders  of  restoration  higher 
than  h = 3 or  4. 

Figure  3 shows  the  phase  error  for  the  operator  (101  as  a function  of  level  of 
restoration  (II  in  the  figure)  and  wavelength  in  units  of  the  original  grid  length  (Ax'. 
The  phase  difference  is  shown  for  r 0.  25.  for  which  the  phase  error  is  close  to  its 
maximum  value  for  ill  wavelengths.  For  h 0,  the  phase  error  is  greater  than 
10  * radians  for  waves  up  to  about  15  to  16  grid  lengths.  For  h - 1.  the  corres- 
ponding error  occurs  for  waves  !!  to  9 gr  id  lengths.  For  h 2.  the  threshold  is 
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reached  at  w ives  of  (i  to  7 grid  lengths  and  for  li  :i,  the  threshold  is  only  about 
a grid  lengths.  Hut  for  ti  - (>,  it  decreases  only  slightly  to  4 grid  lengths.  The 
amplitude  error  for  ( 10)  with  q h is  essentially  the  same  as  that  of  the  amplitude 
restoring  oper  itor  (9)  which  is  illustrated  in  I igure  1 for  r 0.  r>.  Again,  we  see 
[ that,  except  for  short  waves,  little  further  amplitude  restoration  takes  place  for 

h greater  than  :i  or  4. 

: 
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I igure  3.  The  Phase  Krror,  in 
Radians  for  r 0.  25  for  the 
Phase  and  \mplitude  Restoring 
Interpolation  Operator  (10)  as  a 

I unction  of  Level  of  Restoration 

II  and  Wavelength.  The  11=  0 
curve  corresponds  to  the  0.  2a 
curve  in  Figure  2 


For  applications  where  some  damping  of  the  amplitude  of  the  short  waves  is 
tolerable  or  desirable,  there  appears  to  be  no  advantage  in  exceeding  levels  of 
restoration  corresponding  to  h 3 or  4.  It  should  be  noted,  in  view  of  the  fact  that 
the  magnitude  of  both  the  phase  and  amplitude  error  depends  upon  the  position  of  the 
interpolation  point  m between  the  grid  points  i and  is  1,  that  one  mav  wish  to  design 
an  interpolation  procedure  in  which  the  levels  of  restoration  h and  q are  allowed  to 
vary  with  position  so  as  to  compensate  for  this  variation. 

There  is  considerable  flexibility  in  the  manner  in  which  operator  (10)  mav  be 
used,  depending  upon  the  nature  of  the  application.  The  ability  to  correct  for  phase 
shift  and  amplitude  truncation  may  find  a wide  variety  of  applications,  but  the 
application  for  which  the  present  procedure  is  particularly  appropriate  is  one  which 
concerns  the  calculation  of  finite-differences  involving  both  interpolated  and  original 
data  points.  The  phase  and  amplitude  errors  produced  by  interpolation  introduce 
an  additional  error  in  the  finite  differences.  It  may  be  beneficial  to  minimize  such 
errors  by  the  use  of  an  appropriate  high-order  phase  and  amplitude  restoring  inter- 
polation operator. 
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3.  COMI’XKISON  OF  I’ll  \S1  \M(  \MI*l.m  111  ttll  Ill’ll \>l 

rkstokim;  imkkpoi ation 

While  phase  and  amplitude  restoring  interpolation  (10)  is  based  upon  the  simpler 
amplitude  restoring  interpolation  (0),  it  is  at  the  same  time  more  versatile  and 
more  limited  in  scope.  Interpolation  with  (10)  is  defined  onlv  for  those  points 
<*  1 + r between  the  original  grid  points  i and  i + 1 such  that  r takes  on  the  dis- 

crete values  S/2^'+'  where  S = 1,  2,  3,  ...  . (2^**  1 - !).  On  the  other  hand,  inter- 
polation with  (0)  is  defined  for  all  points  r between  i and  i + 1.  f urthermore,  in  (0) 
as  h increases,  the  size  of  the  stencil  increases  uniformly  and  the  amplitude  re- 
sponse monotonicallv  approaches  unity  for  all  wavelengths  except  the  two-grid- 
interval  wavelength.  With  (10),  however,  although  phase  is  increasingly  preserved 
as  h and  q increase  for  stencils  of  the  same  size,  (9)  is  more  efficient  in  amplitude 
preservation.  Because  of  their  differing  characteristics  direct  comparisons  were 
made  between  operators  (9)  and  (10)  in  both  a static  and  a dynamic  framework. 

3.1  Sialic  Comparisons 

[n  the  static  tests  f(x)  is  represented  in  turn  by  a single  periodic  l-'ourier  com- 
ponent varying  from  wave  number  1 through  wave  number  35  in  the  interval  0 to 
27?.  Discrete  values  of  f(x)  are  sampled  at  uniform  subintervals  Ax  = 277/7  2.  These 
discrete  values  serve  as  data  points  from  which  interpolated  values  are  obtained 
within  each  subinterval  at  specific  values  of  r = 0.  125,  0.  25,  0.  375,  and  0.  5.  The 
root  mean  square  error  is  then  obtained  for  each  wave  number  and  for  each  value 
of  r by  comparing  interpolated  values  with  the  functional  values.  The  results  are 
shown  in  Figures  4 through  8. 

In  these  figures,  the  ordinate  RMSE  is  the  root  mean  square  error  for  either 
the  amplitude  or  phase  and  amplitude  restoring  interpolation,  defined  as. 


where  f.  „ is  the  functional  value  of  f(x)  at  the  point  i + R between  i and  i*  1 and  T. 
is  the  corresponding  interpolated  value.  The  RMSE  curves  for  amplitude  restoring 
interpolation  are  indicated  in  the  figures  by  RMSE  followed  by  a single  pair  of 
parentheses.  The  number  in  the  parentheses  indicates  the  number  of  terms  in  (he 
interpolation  stencil.  The  curves  for  the  phase  and  amplitude  restoring  interpola- 
tion are  indicated  by  RMSE  followed  by  two  pairs  of  parentheses.  The  numbers  in 
the  first  pair  of  parentheses  indicate  the  order  of  the  interpolation  operator  (h,  q) 
and  the  number  in  the  second  pair  of  parentheses  indicates  the  number  of  terms  in 
the  interpolation  stencil.  For  each  wave  number  t'(\),  and  therefore  f.^,  is  com- 
pletely determined  by  the  phase  and  amplitude  of  the  wave.  For  each  wave  the  amp! 
tude  was  taken  as  1 and  the  phase  was  taken  as  -kAx/2  where  k is  wave  number. 
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In  Figure  4,  with  r = R 0.  5,  M = 0 and  therefore  the  amplitude  restoring 
interpolation  operator  (9)  is  identical  to  the  phase  and  amplitude  restoring  operator 
(10).  The  RMSE  is  shown  for  f*0Q  5 (the  2-point  operator),  through  fj4^  5 (the 
10 -point  operator).  For  any  wave  number,  RMSE  monotonically  decreases  as  the 
order  of  the  operator  increases.  However,  the  decrease  in  RMSE  is  large  and 
significant  for  the  lower  wave  numbers  and  vanishingly  small  for  the  largest  wave 
numbers.  If  one  considers  the  spectral  energy  density  distribution  of  most  geo- 
physical data,  most  of  the  variance  is  in  the  lower  wave  numbers  and  most  of  the 
noise  is  in  the  highest  wave  numbers.  Therefore,  for  geophysical  applications, 
the  distribution  of  RMSE  with  wave  number  is  highly  advantageous,  especially  for 
the  higher  ordered  operators. 


Figure  4.  Root  Mean  Square  Error 
(RMSE)  for  Amplitude  Restoring 
Interpolation  and  Phase  and  Amplitude 
Restoring  Interpolation  for  Various 
Wave  Numbers  in  a 7 2 -Grid-Interval 
Domain.  The  amplitude  of  each 
wave  is  1.  0 regardless  of  its  wave 
number.  Numbers  in  single 
parentheses  indicate  the  size  of  the 
interpolation  stencil.  When  the 
interpolation  point  R - 0.5,  ampli- 
tude and  phase  and  amplitude 
restoring  interpolation  are  identical 


In  Figure  5,  r * 0.  125,  there  is  a clear  distinction  between  the  amplitude  re- 
storing operators  and  the  phase  and  amplitude  restoring  operators.  The  uppermost 

curve  in  the  figure  corresponds  to  f'0!,  the  2-point  amplitude  restoring 

i‘°.  125  (1,  (2i 

operator.  The  next  curve  below  this  contains  the  superposition  of  E(  o.  12r,-  Vo.  125' 

and  f*4*  it  is  apparent  that  almost  all  of  the  improvement  in  RMSE 

i*0.  125  t*0.  125 

that  can  be  produced  by  amplitude  restoration  alone  is  effected  by  the  4 -point 
operator,  f(  ^ . Further  improvement  in  RMSE  can  be  obtained  bv  correcting 
for  phase  error  at  the  same  time  that  amplitude  restoration  is  taking  place,  that  is. 
bv  using  operator  (10)  rather  than  operator  (!)).  There  is  a steady  decrease  in  the 
size  of  RMSE  with  the  use  of  (10)  as  the  order  of  the  operator  increases  from 
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Figure  5.  Same  as  Figure  4 with 
H 0.  125.  The  order  (h,  q)  of 
phase  and  amplitude  restoring 
interpolation  is  indicated  bv  the 
numbers  in  the  first  pair  of 
parentheses  for  those  curves 
identified  by  two  pairs  of  paren- 
theses. The  number  in  the  second 
pair  of  parentheses  indicates  the 
size  of  the  interpolation  stencil 


h q - 1 through  h q 4.  In  fact,  the  5 -point  operator  obtained  with  the  use  of 
(10),  with  h q 1,  provides  a RMSE  distribution  which  is  superior  to  that  arising 
from  the  use  of  the  10-point  operator  fj^  Similar,  though  not  identical  re- 

sults, are  obtained  at  other  values  of  r as  can  be  seen  from  Figures  6 and  7 which 
show  the  RMSE  distributions  for  r 0.25  and  0.375,  respectively.  In  these  cases, 
RMSE  (1,  1)  is  superior  to  RMSE  (10)  only  for  the  lower  wave  numbers,  but  RMSE 
(2,  2)  is  clearly  superior. 


Figure  6.  Same  as  Figure  5 
with  R - 0.  25 
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In  Figure  8,  for  r 0.25,  a sampling  of  results  is  shown  for  operator  (10)  with 
h 4 q.  Although,  in  general,  the  results  with  h / q are  inferior  to  those  with  h - q, 
even  where  the  stencil  sizes  are  comparable,  these  results  illustrate  one  of  the 
versatile  features  of  the  phase  and  amplitude  restoring  interpolation;  namely,  the 
size  of  the  stencil  can  be  tailored  to  fit  changing  circumstances. 


Figure  8.  Same  as  Figure  5 
with  R = 0.  25 
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3.2  Dynamic  (xitnparisons 

Dynamic  tests  were  carried  out  in  a simple  system  that  approximates  the  com- 
plexity of  operational,  multilevel,  primitive  equation  models.  1 he  eight-level 
primitive -equation  channel  model  is  fully  described  in  an  earlier  study.  It  was 
used  to  compare  amplitude  restoring  interpolation  to  two-point  linear  interpolation, 
with  results  which  are  outlined  in  the  Introduction.  The  same  basic  framework  is 
used  in  the  present  study  with  some  differences  in  detail,  necessitated  primarily  in 
order  to  accommodate  tests  with  phase  and  amplitude  interpolation.  Thus,  in  the 
present  study,  the  fine-mesh  grid  size  is  1/4  the  size  of  the  coarse-mesh  grid 
(555/4  km)  rather  than  1/5  as  in  the  earlier  study.  The  amplitude  restoring  inter- 
polation in  the  present  study  makes  use  of  the  more  efficient  inverse  operator  (6) 

rather  than  the  multiplicative  operator  described  by  the  author.  In  addition,  the 
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smoothing  operators  used  in  the  present  study  are  not  only  more  efficient,  but 
are  applied  in  a slightly  different  fashion.  Finally,  in  order  to  permit  an  unbiased 
basis  for  judging  the  merits  of  each  test  solution,  a fine  mesh  validation  (FMY) 
solution  was  obtained. 

All  test  solutions,  including  FMV,  start  with  the  same  smooth  coarse-mesh 
solution  at  time  step  6000  (83.  33  days)  as  initial  data.  IN  FMV,  the  coarse-mesh 
initial  fields  of  velocity  and  temperature  were  interpolated  throughout  the  entire 
coarse-mesh  domain  using  the  high-order  phase  and  amplitude  restoring  inter- 
polation operator  (10),  with  h q 3 wherever  possible,  and  reducing  the  order  of 

the  operator  only  in  approaching  the  northern  or  southern  boundary.  FMV  solutions  t 

were  obtained  throughout  the  entire  domain  with  the  fine -mesh  solution  in  space 
(138.75  km)  and  time  (5  min),  for  1600  time  steps.  These  solutions  then  serve  as 

validation  solutions  for  the  fine-mesh  nested  solutions  since  they  are  obtained  J 

without  any  artificial  internal  boundaries  and  except  for  the  initial  interpolation  at 
83.33  days,  without  anv  further  interpolation. 

Three  separate  series  of  test  solutions  in  the  nested  fine-mesh  domain  were 
carried  out  from  the  same  coarse-mesh  initial  conditions  and  with  precisely  the 
same  conditions  in  all  respects  except  for  the  type  of  interpolation  used  to  obtain 
the  initial  information  and  continuing  boundary  information. 

Series  A,  amplitude  restoration,  makes  use  of  Fq.  (9)  with  h 3.  Series  B, 
phase  and  amplitude  restoration,  makes  use  of  Fq.  (10),  with  h q 3 at  all  points 
where  it  is  possible  to  do  so.  Since  the  northern  boundary  of  the  fine-mesh  domain 
is  situated  at  17.5  deg  N,  there  are  only  three  coarse-mesh  grid  points  north  of 
this  boundary.  After  the  initialization  for  the  fine-mesh  experiments  interpolation 
is  performed  only  one  dimensionally  from  the  continuing  coarse-mesh  solutions 
onto  the  fine-mesh  boundaries.  Since  the  coarse  mesh  domain  is  periodic  in  the 
east -west  direction,  there  is  no  difficulty  in  using  Fq.  (10),  with  h q 3 on  the 
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northern  and  southern  fine- mesh  boundaries.  However,  it  is  necessary  to  gradually 
decrease  the  order  of  interpolation  on  the  eastern  and  western  boundaries  for  some 
of  the  grid  points  near  the  northern  boundary. 

The  grid  spacing  is  uniform  in  the  east-west  (x)  ana  north-south  directions  and 
is  555  km  in  the  coarse-mesh  and  138.75  km  in  the  tine-mesh  configuration.  The 
vertical  spacing  is  125  mb  in  both  mesh  systems.  The  northern  and  southern 
boundaries  are  located  at  32.  5 deg  N and  32.  5 deg  S in  the  large-scale  model,  and 
at  !7.  5 deg  N and  2.  5 deg  N in  the  fine-mesh  model  (see  Figure  91.  The  east-west 
domain  extends  over  20  grid  intervals  (11,  100  km)  and  is  periodic  in  the  coarse- 
mesh  domain.  In  the  fine-mesh  nested  region  the  east-west  domain  contains  16  grid 
intervals  (2220  km)  and  is  not  periodic.  The  north-south  domain  contains  13  grid 
intervals  in  the  coarse-mesh  region  (7215  km'  and  12  grid  intervals  (1665  km'  in 
the  fine-mesh  nested  region.  Thus,  if  i and  j represent  east-to-west  and  south-to- 
north  grid  point  index  designations,  i varies  from  1 through  17  and  j,  from  1 through 
13  in  the  fine-mesh  nested  region.  Thus,  at  j = 8,  north-south  interpolation  on  the 
eastern  and  western  boundaries  makes  use  of  a lower-order  form  of  Eq.  (10)  with 
h = 2 and  q = 3.  At  j = 0,  on  the  eastern  and  western  fine-mesh  boundaries  there  is 
no  need  for  interpolation  since  these  points  correspond  to  coarje-mesh  points.  This 
is  also  true  at  j = 13.  At  j = 10  we  use  Eq.  (10)  with  h = 3 and  q = 2 and  at  j = 12, 
with  h = q = 2.  At  i = 11,  the  interpolation  point  a lies  midway  between  the  coarse- 
mesh  points  and  so  Eq.  (10),  with  h q 3,  reduced  to  Fq.  (0),  with  h 3. 

In  the  Series  C tests,  only  two-point  linear  interpolation  is  used. 

To  recapitulate,  the  following  procedure  was  followed  in  the  various  tests. 

(1)  Starting  from  a stationary  motion  field,  with  no  horizontal  temperature 
differences  and  a vertical  temperature  distribution  appropriate  for  the  tropics,  the 
heating  function  is  turned  on  in  the  large-scale  domain  and  motion  begins  to  develop. 

(2)  The  coarse-mesh  run  is  continued  beyond  6400  time  steps  of  20-min  duration 
for  a total  elapsed  time  of  about  90  days. 

(3)  Starting  at  time  step  3000  we  performed  occasional  smoothing  with  a high- 
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order,  amplitude  restoring  operator  on  the  predicted  fields  of  horizontal  velocity 
and  temperature.  The  9-point  smoothing  operator,  shown  in  Table  2 with  n 3,  was 
used  in  both  the  east-west  and  north-south  directions  whenever  it  was  possible  to  do 
so.  For  those  points  situated  four  grid  points  or  less  from  the  northern  or  southern 
boundary  7-,  5-,  and  3-point  operators  were  used  as  appropriate  in  both  the  east- 
west  and  north-south  directions.  These  operators  are  also  shown  in  Table  2 with 
p = 2,  1,  or  0,  respectively.  Thus,  although  the  smoothing  was  performed  only 
occasionally,  the  boundary  regions  of  the  domain  w ere  more  strongly  smoothed 
inasmuch  as  lower-order  operators  were  used  near  the  boundaries.  The  smoothing 
was  performed  only  every  1000  time  steps  between  time  steps  3000  and  5000,  but 
then  the  frequency  of  smoothing  was  increased  to  every  10  time  steps  from  step  5010 
through  step  6000.  This  produced  smoothed  fields  for  the  initial  times  for  the  various 
test  runs. 
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east;  'lashed  curves  from  east  to  west.  The  rectangular  \ns 
of  the  fine-mesh,  nested  domain 


for  Various  Values  of  p 


Table  2.  Stencils  for  f.*^" 


p 

f. 

1 

r ■ 

fi*l 

fi±  2 

fi±  3 

fi±4 

fi*6 

fi±7 

fi±  8 

fi±  9 

fi±  10 

0 

1/22 

(2 

1) 

1 

1/24 * 

(10 

4 

-1) 

2 

1/26 * 

(44 

15 

-6 

1) 

3 

1/28 

(186 

56 

-28 

8 

-1) 

4 

1 / 2 1 0 

(772 

210 

- 120 

45 

- 10 

1) 

5 

1 /2 1 2 

(3172 

7 92 

-495 

220 

-66 

12 

-1) 

6 

1 /2 14 

(12952 

3003 

-2002 

1001 

-364 

91 

-14 

1) 

7 

1 / 2 1 6 

(52666 

1 1440 

-8008 

4368 

- 1820 

560 

- 120 

16 

-1) 

8 

1/218 

(213524 

43758 

-31824 

18564 

-8568 

-816 

153 

-18 

1) 

9 

1/220 



(863820 

167960 

— 

- 125970 

-4845 

1 140 

- 190 

20) 

-1) 

(4)  For  the  principal  series  of  experiments,  starting  from  day  83.  33,  both  the 

coarse-mesh  solutions  (which  were  used  to  supply  boundary  information  for  the 

nested,  fine-mesh  solutions)  and  the  various  fine-mesh  solutions  (including  FMV) 
were  smoothed  at  2000-min  intervals  (every  100  coarse-mesh  time  steps  or  every 
400  fine-mesh  time  steps).  This  is  a departure  from  the  earlier  study'*'1  in  which 
the  smoothing  interval  was  only  200  min.  However,  other  tests  in  the  present  study 
were  carried  out  with  more  frequent  smoothing;  namely,  at  intervals  of  250  min  and 

50  min  in  order  to  separate  the  effects  of  smoothing  from  those  of  i nterpolation. 

The  results  of  the  various  test  solutions  and  comparisons  with  the  FMV  solutions 
are  shown  for  two  separate  times;  namely,  at  800  fine-mesh  time  steps  (day  86.  11) 
and  at  fine-mesh  time  step  1600  (day  88.  89).  As  in  the  earlier  study,  the  east- 
west  component  of  the  velocity  (u)  is  used  to  illustrate  the  relative  merits  of  the 
various  test  runs.  The  results  are  shown  both  graphically  in  several  figures  and 
in  Tables  3 and  4,  which  show  the  correlations  between  FMV  and  the  various  test 
solutions  as  well  as  the  root  mean  square  differences  between  them. 
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Table  3.  Root  Mean  Square  Differences  (RMSD)  of  the  Zonal  Wind  Field  in  cm  sec 
Between  Fine  Mesh  Validation  Smoothed  Every  2000  min  (FW  20001  and  other 
Model  Solutions  After  800,  5-minTime  Steps  (day  86.  Ill  and  1600  Time  Steps 
(day  88.  89).  Results  apply  only  to  the  limited-area  domain 


Level 

Day  86 

. 1 1 

AS2000 

BS2000 

CS2000 

AS250 

B.S250 

CS250 

\S50 



BS50 

CS50 

FMV250 

n 

25.  3 

24.  4 

31.7 

9.4 

9.  3 

12.  8 

5.  5 

5.  5 

9.  2 

0.  3 

H 

12.  8 

12.4 

13.6 

8.  9 

8.  9 

10.  3 

8.  5 

8.  5 

10.  2 

0.  4 

I 

7.  1 

6.  8 

7.  1 

5.  5 

5.  5 

6.  0 

5.  8 

5.  8 

6.  6 
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Figure  10  shows  the  distribution  of  zonal  velocity  in  FMV  at  day  86.  11  for  two 
upper  levels  (level  6,  312.5  mb;  level  8,  437.5  mb)  and  two  lower  levels  (level  12, 
687. 5 mb;  level  14,  812.  5 mb),  but  only  for  the  area  occupied  by  the  nested  region. 
Figures  11  through  13  show  the  comparable  results  for  Series  A,  B,  and  C.  Flows 
from  west  to  east  are  indicated  by  solid  lines  and  flows  from  east  to  west  by  dashed 
lines.  In  Figures  10  through  13,  smoothing  is  performed  every  400  time  steps, 
that  is,  only  twice  from  the  initial  time  (day  83.  33)  until  day  86.  11.  In  spite  of  the 
insignificant  amount  of  smoothing  the  FMV  results,  for  the  nested  region,  are 
smooth  at  all  levels.  However,  the  results  shown  for  Series  A,  B,  and  C are  not 
as  consistently  smooth.  There  is  little  difference  among  the  three  series  at  the 
two  upper  levels,  and  each  shows  quite  good  agreement  with  FMV,  at  least  as  far 
as  the  larger  scale  features  are  concerned.  The  minor  differences  that  are  appar- 
ent among  the  three  series  are  in  the  direction,  not  unexpectedly,  of  increasing 
smoothness  and  versimilitude  (with  regard  to  FMV)  in  proceeding  from  Series  O 
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Table  4.  I. inear  Correlation  Coefficients  Between  the  FMV  2000  Zonal  \\  ind  Field 
and  the  Zonal  Wind  Fields  for  Other  Model  Solutions  Calculated  for  the  Limited- 
Area  Domain 


Level 

Day  8 b 

. 1 1 

AS2000 

BS2000 

CS2000 

AS  250 

BS250 

CS250 

AS50 

BS50 

CS50 

FMV  250 

2 

. 991 

. 992 

. 986 

. 999 

. 999 

. 999 

1.  000 

1.  000 

1.  000 

1.  000 

4 

. 990 

. 992 

. 990 

. 998 

. 998 

. 999 

. 998 

. 998 

. 999 

1. 000 

6 

. 944 

. 949 

. 937 

. 969 

. 969 

. 957 

. 959 

. 959 

. 945 

. 999 

8 

. 946 

. 951 

. 947 

. 981 

. 982 

. 97  9 

. 986 

. 985 

. 980 

. 999 

10 

. 889 

. 891 

. 885 

.912 

. 909 

. 929 

. 928 

. 929 

. 938 

. 995 

12 

. 952 

. 95  l 

. 954 

. 970 

. 96  8 

. 975 

. 97 6 

. 97  5 

. 978 

. 999 

14 

. 983 

. 984 

. 97  9 

. 994 

. 994 

. 994 

. 995 

. 995 

. 995 

1.  000 

16 

. 974 

. 975 

. 974 

. 984 

. 984 

. 987 

. 986 

. 985 

. 988 

1.  000 

Day  88.  89 

2 

. 97  8 

. 979 

. 97  0 

. 998 

. 998 

. 998 

. 999 

. 999 

. 999 

1.  000 

4 

. 982 

. 984 

. 983 

. 996 

. 996 

. 997 

. 996 

. 996 

. 996 

1.  000 

6 

. 685 

. 702 

. 704 

. 771 

. 772 

. 773 

.783 

.785 

o 

00 

I- 

. 995 

8 

. 866 

. 869 

. 874 

. 919 

. 915 

. 928 

. 933 

. 932 

. 927 

. 996 

10 

. 557 

. 554 

. 574 

. 631 

. 631 

. 647 

. 623 

. 621 

.632 

. 994 

12 

.810 

. 817 

. 773 

. 869 

. 869 

. 860 

. 866 

. 866 

. 865 

. 998 

14 

. 968 

. 96  9 

. 961 

. 984 

. 983 

. 987 

. 984 

. 985 

. 987 

1. 000 

16 

. 965 

. 966 

. 958 

. 988 



. 988 

. 991 

. 989 

. 989 


. 991 


1.  000 

to  Series  A to  Series  B.  At  the  lower  levels,  however,  there  are  sizable  and  evident 
boundary  effects  in  the  C series  two-point  interpolation',  but  relatively  minor 
boundary  effects  in  the  A and  B series.  On  the  other  hand,  the  results  from  the  A 
and  B series  are  virtually  identical.  In  spite  of  the  obvious  boundary  effects  in 
Series  C at  the  lower  levels,  the  larger  scale  features  of  all  three  series  closely 
resemble  the  FMV  results. 

The  above  conclusion  is  also  evide-*  from  Tables  3 and  4 showing  the  root  mean 
square  differences  (RMSD)  and  linear  correlation  coefficients,  level  for  level,  be- 
tween FMV  and  each  of  the  three  series.  In  Tables  3 and  4,  the  series  of  solutions 
(A,  B,  C)  as  well  as  the  time  intervals  between  smoothing  are  indicated  bv  the 
column  headings.  Thus,  AS2000  indicates  Series  A with  smoothing  every  2000  min- 
utes. 


Figure  11.  Same  as  Figure  10  for  Series  A,  Amplitude  Restoring  Interpolation 


It  is  apparent  that  there  is  little  difference  among  the  solutions  AS2000,  BS2000, 
and  CS2000  both  with  regard  to  the  root  mean  square  differences  as  well  as  with 
regard  to  the  correlation  coefficients.  It  is  also  apparent,  as  far  as  the  RMSD's 
and  correlation  coefficients  are  concerned,  that  it  is  the  large-scale  features  that 

contain  most  of  the  variance  and  therefore  dominate  the  results. 
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In  the  earlier  study,  ' where  1/5  mesh  size  was  used  in  experiments  compar- 
able to  Series  A and  Series  C,  there  were  far  more  pronounced  boundary  effects 
especially  where  two-point  interpolation  was  used.  The  more  evident  boundary 
effects  occurred  in  the  earlier  study  in  spite  of  much  more  frequent  smoothing 
(every  200  min  as  compared  with  every  2000  min  in  the  present  study).  This  re- 
sult points  up  the  prominent  role  played  by  phase  error  in  the  boundary  information 
since  with  1/5  mesh  size  there  is  appreciable  phase  error  at  all  interpolated  points, 
whereas  with  1/4  mesh  size  there  is  no  phase  error  at  the  central  interpolated  point 
even  in  Series  A and  Series  C.  Nevertheless,  it  is  apparent  from  f igure  12  that 
correcting  for  both  phase  and  amplitude  error  in  the  interpolated  boundary  informa- 
tion is  not  sufficient  by  itself  to  avoid  boundary  effects.  This  conclusion  is  even 
more  apparent  from  the  results  of  the  solution  after  1600  time  steps  (day  88.  8P>. 

Figures  14  through  17  are  comparable  to  Figures  10  through  13  but  illustrate  the 
results  for  day  88.  89.  The  Series  A and  B results  are  somewhat  noisier  than  they 
were  on  day  86.  11,  particularly  at  the  upper  two  levels.  Series  C,  Figure  17,  how- 
ever, is  significantly  noisier,  at  all  levels  on  day  88.  89  as  compared  with  the  result  s 
on  day  86.  11.  However,  in  all  three  series  the  large-scale  distribution  has  not 
changed  appreciably.  The  F MV  results  on  day  88.  89  remain  smooth,  but  show  a 
new  smaller-scale  variation  in  the  north-south  direction  which  was  not  evident  on 
day  86.  11.  None  of  the  three  test  series  (A,  B,  or  C)  show  any  evidence  of  this 
smaller-scale  variation,  although  the  larger-scale  features  of  each  continue  to 
resemble  the  larger-scale  features  of  FMY  on  day  88.  89.  This  continued  resem- 
blance of  the  larger-scale  features  is  reflected  in  Tables  3 and  4 by  the  still  rela- 
tively small  RMSD's  and  relatively  large  correlations  between  each  of  the  series  and 
FMV  on  day  88.  89.  There  is  some  growth  of  RMSD  and  some  small  decrease  in 
correlation  between  days  86.  11  and  88.  89,  especially  at  the  mid-levels  where  the 
flow  changes  from  predominantly  easterly  to  predominantly  westerly.  However,  in 
spite  of  the  growth  of  noise  in  the  three  test  series,  it  is  apparent  that  the  RMSD's 
and  correlations  are  dominated  by  the  variance  contained  in  the  large-scale  features. 

I shall  return  subsequently  to  the  question  of  the  smaller-scale  features  in  1M\  on 
day  88.  89  which  are  not  captured  by  any  of  the  test  series. 
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It  is  apparent  from  the  results  that  minimizing  phase  error  in  the  bound  a r\ 
information  used  to  obtain  the  limited  area  solutions  is  an  important  aspect  in  main- 
taining reasonable  agreement  between  the  test  solution  and  the  validation  solution. 
However,  it  is  also  apparent,  is  has  been  mentioned,  that  there  is  n increasing 
growth  of  boundary  generated  noise  as  the  solutions  are  carried  forward  in  time, 
particularly  in  the  Series  ('  solutions.  In  the  earlier  study,  with,  fine-mesh  spacing 
1 “>  the  coarse-mesh  spacing,  smoothing  of  the  fields  was  carried  out  much  more 
frequently  (every  200  mini  and  vet,  the  results  were  affected  to  n eve1:  greater 
extent  by  boundary  generated  noise.  In  order  to  determine  whether  more  frequent 
smoothing  of  the  present  solutions  would  be  effective  in  lecre  sing  tin-'  noise,  two 
additional  sets  of  solutions  were  obtained  with  the  Series  A,  B,  nd  i conditions. 

In  each  case  the  new  Series  A,  B,  and  C solutions  were  obtained  in  precisely  the 
same  manner  as  the  original  series  except  for  the  frequency  of  sr:  oothing.  In  one 
set  the  smoothing  was  performed  every  250  min  and  in  the  other,  every  50  mm. 

The  IlMSD's  and  correlations  between  these  -ets  of  solutions  old  the  original  I \1\ 
solution  are  tabulated  in  Tables  3 and  4,  under  columns  1.  beled  ''2  40  or  40.  I lu 
zonal  wind  fields  at  levels  6,  8,  12,  and  14  are  shown  for  the  AS240,  I . '250,  uid 
CS250  for  day  86.  11  (Figures  18  through  20)  and  day  88.8°  il  igures  21  through  22). 
Figures  24  through  26  and  Figures  27  through  20  show  the  comp  ir  hie  results  for 
Series  A,  B,  and  C'  with  smoothing  every  50  min  (that  is,  every  ten  time  steps). 

It  is  immediately  apparent  from  the  figures  that  the  additional  smoothing,  even 
at  intervals  of  250  min,  is  sufficient  to  smooth  out  the  minor  boundary  induced  ir- 
regularities in  all  three  series  at  the  two  upper  levels,  f urthermore,  the  major 
irregularities  in  the  C series  are  also  smoothed  out  at  the  two  lower  levels.  In 
addition,  the  increased  frequency  of  smoothing  has  the  effect  of  substantially  reduc- 
ing the  KMSD's  and  increasing  the  linear  correlation  between  all  three  series  and 
1 MV.  Compared  with  the  earlier  study  in  which  smoothing  was  performed  every 
200  min,  the  improvement  in  all  three  series,  with  smoothing  every  250  min,  in 
reducing  boundary  induced  irregularities,  decreasing  the  size  of  IlMSD's  and  in- 
creasing the  correlation  between  the  series  and  FMV  is  indeed  remarkable  and 
igain  focuses  attention  on  the  important  difference  between  a 1 5 mesh  size  and  the 
1 4 mesh  size  of  the  present  study.  The  difference  between  the  former  1 5 and  the 
present  1/4  mesh  sizes  is  the  only  important  difference  between  the  earlier  and  the 
present  studies  as  far  as  Series  AS250  and  CS250  (amplitude  restoring  ,nd  two-point 
interpolation!  are  concerned. 


49 


Smoothing  Every  250 


Figure  21.  Same  as  Figure  15  with  Smoothing  livery  250  nun 


Smoothing 


Same  as  Figure  12  with  Smoothing  I very  50  min  (10  tune  steps) 


Smoothing  Every  50  min  (10  time  steps) 


Smoothing 


29.  Same  as  Figure  17  with  Smoothing  Every  50  min  (10  time  steps) 


Judging  from  the  figures  as  well  as  from  Tables  3 and  4,  all  three  series 
generally  show  some  improvement  with  smoothing  every  50  min  as  compared  with 
those  obtained  with  smoothing  every  250  min.  However,  the  differences  are  rather 
insignificant  except  in  the  case  of  the  RMSD's  at  level  2.  The  results  for  the  B 
Series  are  generally  superior  to  those  of  the  A and  C Series,  but  here  too,  the 
differences  are  small.  It  is  conceivable,  with  different  initial  conditions  especially 
where  a greater  fraction  of  the  variance  is  contained  in  smaller-scale  features,  that 
there  might  have  been  a larger  difference  between  the  results  with  smoothing  every 
250  min  and  every  50  min  or  among  the  results  for  Series  A,  B,  and  C.  However, 
the  results  of  the  present  study  in  conjunction  with  those  of  the  earlier  study  show 

(1)  that  it  is  of  vital  importance  to  control  and  minimize  the  phase  error  that 
arises  from  interpolation  in  obtaining  boundary  information  in  a fine-mesh,  limited- 
area  solution; 

(2)  that  if  the  phase  error  is  controlled,  a relatively  small  amount  of  smoothing 
is  sufficient  to  yield  reasonably  smooth  and  accurate  results; 

(3)  that  the  amount  of  smoothing  necessary  to  yield  smooth  and  accurate  re- 
sults depends  directly  on  the  care  taken  to  control  phase  error  on  the  boundaries; 

(4)  that,  therefore,  it  is  preferable  to  control  phase  error  on  the  boundaries 
by  using  a fine-mesh  grid  size  that  maintains  a ratio  of  (l/2)n  of  the  coarse-mesh 
grid  in  conjunction  with  a relatively  simple  interpolation  scheme  such  as  was  used 
in  Series  A. 

4.  ADDITIONAL  RESULTS 

It  had  been  noted  that  by  day  88.  89,  the  FMV  solution  develops  a small-scale 
variation  that  does  not  develop  in  any  of  the  test  solutions  (Series  A,  B,  or  C)  or 
for  that  matter,  in  the  coarse-mesh  solution  that  is  used  to  supply  the  boundary 
information  for  the  test  solutions.  Although  the  FMV  solution  remains  smooth  on 
day  88.  89,  it  is  possible  that  the  small-scale  feature  is  a numerically  generated 
artifact  rather  than  the  result  of  a physical  process  that  becomes  manifest  in  FMV 
only  because  of  the  increased  resolution.  To  help  resolve  this  uncertainty  an  addi- 
tional experiment  was  performed.  A new  FMV  solution  was  obtained  with  smooth- 
ing every  250  min.  A comparison  between  the  new  FMV  250  and  the  original  FMV 
2000  is  shown  in  Tables  3 and  4 under  the  column  heading  FMV  250.  It  is  apparent 
from  the  small  HMSD's  and  the  large  correlation  coefficients  that  the  solutions  are 
virtually  identical.  It  is  also  apparent  from  Figures  30  and  31,  which  show  the 
FMV  250  results  that  are  comparable  to  the  FMV  2000  results  of  Figures  10  and 
14,  that  there  is  scarcely  any  noticeable  difference  between  the  solutions.  It  thus 
seems  highly  unlikely  that  the  small-scale  variation  in  FMV  on  day  88.  89  is  a 
numerical  artifact,  since  such  artifacts  typically  develop  first  in  the  shortest  wave- 
lengths which  would  have  been  severely  dumped  by  the  more  frequent  smoothing. 
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Figure  30.  Same  as  Figure  10  with  Smoothing  Every  250  min  (50  time  steps) 


5.  CONCI.I  SIGNS 


The  simple,  pragmatic  procedure  proposed  here,  namely,  overspecification 
of  the  artificial  lateral  boundaries  of  the  limited-area  domain,  in  conjunction  with 
(1)  care  to  ensure  that  approximately  correct  phase  information  is  supplied  to  the 
boundaries,  and  (2)  a weak  but  high-order  smoothing  operation,  appears  to  yield 
satisfactory  fine-mesh  limited-area  solutions  in  a nested  region  both  in  terms  of 
controlling  noise  and  affording  verisimilitude  to  a fine-mesh  validation  solution. 
Because  of  the  simplicity  of  the  procedure  it  appears  to  have  applicability  to 
operational,  fine-mesh,  limited-area  models.  However,  the  results  also  indicate 
that  unless  more  information  is  contained  in  the  initial  data  of  the  limited  area  than 
is  available  solely  from  the  large-scale,  coarse-mesh  model,  the  fine-mesh  solu- 
tion is  not  likely  to  differ  appreciably  from  what  would  have  been  obtained  from  the 
coarse  mesh  alone. 
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